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The Einstein@Home search for periodic gravitational waves in LIGO S4 data 
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A search for periodic gravitational waves, from sources such as isolated rapidly-spinning neutron 
stars, was carried out using 510 hours of data from the fourth LIGO science run (S4). The search 
was for quasi-monochromatic waves in the frequency range from 50 Hz to 1500 Hz, with a linear 
frequency drift / (measured at the solar system barycenter) in the range —f/r < f < 0.1 //r, 
where the minimum spin-down age r was 1000 years for signals below 300 Hz and 10 000 years above 
300 Hz. The main computational work of the search was distributed over approximately 100 000 
computers volunteered by the general public. This large computing power allowed the use of a 
relatively long coherent integration time of 30 hours, despite the large parameter space searched. 
No statistically significant signals were found. The sensitivity of the search is estimated, along with 
the fraction of parameter space that was vetoed because of contamination by instrumental artifacts. 
In the 100 Hz to 200 Hz band, more than 90% of sources with dimensionless gravitational wave strain 
amplitude greater than 10~^^ would have been detected. 

PACS numbers: 04.80.Nn, 95.55.Ym, 97.60.Gb, 07.05.Kf 



I. INTRODUCTION 

Gravitational waves are a fundamental prediction of 
Einstein's General Theory of Relativity [1, 2 . But these 
waves are very weak, so although there is compelling in- 
direct evidence for their existence [3 , direct detection has 
so far not been possible. 

In the past decade, advances in lasers, optics and con- 
trol systems have enabled construction of a new gener- 
ation of gravitational-wave detectors [4] that offer the 
first realistic promise of a direct detection. The Laser 
Interferometer Gravitational-wave Observatory (LIGO) 
[5j[6| is currently the most sensitive of these instruments. 
LIGO consists of three kilometer-scale instruments. Two 
are located in a common vacuum envelope in Hanford, 
Washington, USA and the other is located in Livingston, 
Louisiana, USA. 

This paper reports on the results of the Ein- 
stein@Home search for "continuous wave" sources in the 
data from the fourth LIGO science run (S4). The con- 
figuration of the LIGO detectors during the S4 run is 
described in a separate instrumental paper [T. 



Continuous Avave (CW) sources and detection 
methods 



"Continuous waves" (CW) are quasi- monochromatic 
gravitational-wave signals whose duration is longer than 
the observation time. They have a well-defined frequency 
on short time-scales, which can vary slowly over longer 
times. These types of waves are expected, for example, 
from spinning neutron stars with non-axisymmetric de- 
formations. If the system is isolated, then it loses angular 
momentum to the radiation. The spinning slows down, 
and the gravitational- wave frequency decreases. Gravi- 
tational acceleration towards a large nearby mass distri- 
bution can also produce such a frequency drift (of either 
sign). Many possible emission mechanisms could lead 
the to the emission of such waves by spinning neutron 



stars [a[9l|T0l[IIl[ia[T3l[Tl|T5l|T6]. 

If there were no acceleration between the LIGO de- 
tectors and the GW sources, then it would be possible 
to search for CW signals from unknown sources using 
only "standard" computing resources, such as a high- 
end workstation or a small computing cluster. In this 
case the analysis technique would be simple: compute 
the Fast Fourier Transform (FFT) [17, 18 of the original 
time-series data, and search along the frequency axis for 
peaks in the power spectrum. Time-domain resampling 
or similar techniques could be used to compensate for the 
effects of a linear-in-time frequency drift. 

However, this simple analysis is not possible because 
of the terrestrial location of the LIGO detectors: sig- 
nals that are purely sinusoidal at the source are Doppler- 
modulated by the Earth's motion and thus are no longer 
sinusoidal at the detector. The Earth's rotation about 
its axis modulates the signal frequency at the detector 
by approximately one part in 10^, with a period of one 
sidereal day. In addition, the Earth's orbit about the 
Sun modulates the signal frequency at the detector by 
approximately one part in 10^, with a period of one year. 
These two modulations, whose exact form depends upon 
the precise sky location of the source, greatly complicate 
the data analysis when searching for unknown sources. 
The search becomes even more complicated if the CW- 
emitter is part of a binary star system, since the orbital 
motion of the binary system introduces additional mod- 
ulations into the waveform. 

The "brute force" approach to the data analysis prob- 
lem would employ matched filtering, convolving all avail- 
able data with a family of template waveforms corre- 
sponding to all possible putative sources. The resulting 
search statistic is called the JT-statistic and was first de- 
scribed in a seminal paper of Jaranowski, Krolak, and 
Schutz [19 . But even for isolated neutron stars (i.e. 
which are not in binary systems) the parameter space 
of possible sources is four-dimensional, with two param- 
eters required to describe the source sky position using 
standard astronomical equatorial coordinates a (right as- 



cension) and S (declination), and additional coordinates 
(/, /) denoting the intrinsic frequency and frequency 
drift. To achieve the maximum possible sensitivity, the 
template waveforms must match the source waveforms 
to within a fraction of a cycle over the entire observation 
time (with current detectors this is months or years). So 
one must choose a very closely spaced grid of templates in 
this four-dimensional parameter space, and the computa- 
tional cost exceeds all available computing resources on 
the planet [20]. Thus the direct approach is not possible 
in practice. 

More efficient and sensitive methods for this type of 
search have been studied for more than a decade and are 
under development [211 El [23l [Ml ES] . In this paper, the 
frequency-domain method described in [26l [12 is used to 
calculate the ^-statistic. In order to maximize the possi- 
ble integration time, and hence achieve a more sensitive 
coherent search, the computation was distributed among 
approximately 10^ computers belonging to ~ 5 x 10^ vol- 
unteers in ~ 200 countries. This distributed computation 
project, called Einstein@Home [28 , follows the model 
of a number of other well-known volunteer distributed 
computing projects such as SETI@home [29 and Fold- 
ing@home [30 . 

Other methods have also been employed for the CW 
search of the S4 data |3ll [32] and searches for other signal 
types (burst, inspiral, stochastic background) have also 
been carried out [33l [Ml IMl IMIISH [38] with this data 
set. The results of these searches are all upper bounds, 
with no detections reported. 
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FIG. 1: Strain amplitude spectral densities \^Sh{f) of the 
S4 data from the LIGO detectors HI (top) and LI (bottom). 
The gray curves are medians of the entire uncleaned LIGO S4 
science- mode data set with a frequency resolution of 0.125 Hz. 
The black curves show the cleaned S4 data used in this anal- 
ysis with a frequency resolution of 0.5 Hz. The top (bottom) 
plot is the mean of the 10 HI (7 LI) 30-hour data segments 
used in this Einstein@Home analysis. 



B. Outline of this paper 



The outline of this paper is as follows. Sections [ll| 



and III describe the overall construction of the search, 



including the data set preparation, regions of parame- 
ter space searched, and the choices of thresholds and 
sensitivities. Section IV describes the post-processing 
pipeline. The level of sensitivity of the search is esti- 



mated in Section |V[ Section VI describes the vetoing of 
instrumental line artifacts and the fraction of parame- 



ter space that was therefore excluded. Section |VII| de- 
scribes the end-to-end validation of the search and the 
post-processing pipeline, which was done by injecting 
simulated CW signals into the detector hardware. Sec- 
tion [VTTI| describes the final results of the search, followed 
by a short conclusion. 



II. DATA SELECTION AND PREPARATION 

The data for the S4 run was collected between Febru- 
ary 22, 2005 and March 23, 2005. The data analyzed 
consisted of 300 hours of data from the LIGO Hanford 
4-km (HI) detector and 210 hours of data from the LIGO 
Livingston 4-km (LI) detector. 



The search method used here (explained in detail in 
Section III) consists of computing a coherent ^-statistic 



over data segments of 30 hours each, and combining these 
results via an incoherent coincidence scheme. However, 
the 30-hour segments have time-gaps, and the number of 
templates needed for the coherent ^-statistic step grows 
rapidly as the gaps get longer. For this reason, the start 
and end times of the data segments were selected based 
on the criteria that the gaps totaled no more than ten 
hours: each data segment contains 30 hours of science- 
mode data and lies within a total time span of less than 
40 hours. Here and in the following the term "segment" is 
always used to refer to one of these time stretches of data, 
each of which contains exactly Tobs = 30 h of data. The 
total time spanned by the data in segment j is written 

Tspan,,;30h<Tspan,, < 40 h. 

The data segments consist of uninterrupted blocks of 
1800 s of contiguous science-mode data. This is for tech- 
nical reasons: the ^-statistic code uses Short Fourier 
Transforms (SFTs) over Tsft = 1800 s as input data, 
(this data format is described in [39 ). To produce these 
SFTs, the data is first calibrated in the time domain us- 
ing the method described in [40l [41] . Then the data is 
windowed in 1800 s intervals using a Tukey window with 
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TABLE L Segments of S4 data used in this search, in order 
of decreasing sensitivity at 141.3 Hz for HI and at 135.3 Hz 
for LI. The columns are the data segment index j, the GPS 

^^^ and the time spanned 



start time tj, the GPS end time t'j 



a characteristic turn-on/ turn-off time of 500 ms, followed 
by an FFT. 

Applying the above constraints to the S4 data set 
yielded a total of A^seg = 17 data segments (10 from HI, 7 
from LI), labeled by j = 1, • • • ,17. The GPS start time 
of segment j is denoted tj , and these values are listed in 
Table [D 

The maximum Doppler modulation (from the Earth's 
motion about the Sun) is about one part in 10^. Over the 
length of S4, and in the parameter range considered, the 
frequency changes due to intrinsic spin-down are smaller 
still. This means that the CW signals searched for here 
always stays within a narrow frequency band, drifting no 
more than about ±0.15 Hz from some fiducial frequency. 
For this reason the input data, spanning the frequency 
range of 50 Hz to 1500 Hz, is partitioned in the frequency 
domain into 5800 "slices" of 0.5 Hz plus wings of 0.175 Hz 
on either side. The size of one such input data slice is 
7368 000 bytes for HI (containing 600 SFTs from 10 seg- 
ments) and 5 157600 bytes for LI (containing 420 SFTs 
from 7 segments). 

The detector data contains dozens of narrow-band 
spectral lines whose origin is instrumental, for example 
the harmonics of the 60 Hz mains frequency, and vio- 
lin modes of the mirror suspensions in the range from 
342 Hz - 350 Hz (HI) and 335 Hz - 355 Hz (LI). To sim- 
plify later analysis, line features that are known to be 
instrumental artifacts are removed ("cleaned") from the 
data by replacing the frequency-domain data bins with 
computer- generated random Gaussian values. The fre- 
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TABLE H: Instrumental lines removed from the input data. 
The three columns show the frequency of the fundamental 
harmonic /Line, the number of harmonics N, and the band- 
width A/Line removed on either side of the central frequency 
(total bandwidth removed per harmonic is 2 A/Line)- In to- 
tal 77.92 Hz of HI data and 144.29 Hz of LI data have been 
excluded ab initio. If A/Line = then the line-cleaning algo- 
rithm replaces a single Fourier bin with the average of bins on 
either side. The spacing between Fourier bins is 1/1800 Hz. 



quencies of these lines are shown in Table |TI| The clean- 
ing algorithm uses a moving-in-frequency median of the 
power in individual frequency bins to determine the in- 
strumental noise floor. To prevent bias at the boundaries 
of the cleaned regions, the mean of the random values to 
replace the line features interpolates linearly between the 
noise floor at either side of the line feature. The median 
noise strain amplitude spectra of the final cleaned HI and 
LI data sets are shown in Figure [l] 



III. 



DATA PROCESSING 



Figure |2] is a schematic flow-diagram of the Ein- 
stein@Home data processing which is described in this 
section and in the following section on post-processing. 
It shows what parts of the analysis were done by project 
participants, what parts were done on project servers, 
and the relationships between these. 



A. BOINC workunit distribution and validation 

The computational work of the search is partitioned 
into 6 731 410 workunits (separate computing tasks) and 
processed using the Berkeley Open Infrastructure for 
Network Computing (BOINC) [4^, "IS", "IT . Because the 
work is done on computers that are not owned or con- 
trolled by our scientific collaboration or institutions, any 
individual result could be wrong. Error sources include 
defective hardware (such as over-clocked memory) , defec- 
tive software (erroneous system libraries), and malicious 
users (faking correct results). To identify and eliminate 
such errors, BOINC was configured so that each worku- 
nit is done independently by computers owned by at least 
three different volunteers. 

The most common types of errors (lack of disk space, 
corrupted or missing input files, inconsistent internal 
state, etc.) are detected during program execution. If an 
error is detected during run-time, the program reports 
to the Einstein@Home server that the workunit was un- 
successful^ and BOINC generates another instance of the 
workunit, to be sent to another volunteer's computer. 
This behavior is repeated as necessary until three suc- 
cessful results have been obtained. 

The three successful results obtained for each workunit 
are then compared by an automatic validator, which re- 
jects results that do not agree closely. The validation pro- 
cess is more complicated than simple byte-by-byte com- 
parison of output files, because Einstein@Home supports 
multiple computing platforms (Windows, GNU/Linux, 
Mac OS X on Intel and PPC, FreeBSD, and Solaris) and 
differences in CPU hardware, compiler instruction order- 
ing, and floating-point libraries mean that correct and 
valid result files may exhibit small numerical differences. 
The automatic validation takes place in two steps. 

The output files have a fixed five-column format and 
contain 13 000 candidate events, with one line per candi- 
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FIG. 2: Schematic overview of the Einstein@Home data- 
processing and subsequent post-processing. 



date event, as described in Section jIII CI Tfie first vali 



dation step checks that the file syntax is correct and that 
each value is within the allowed range for that column. 
This detects most file corruption. 

Then the validator does comparison of all possible pairs 
of result files. For a given pair of result files, the valida- 
tor checks that corresponding candidate events lie on the 



same template grid-point and have JT-statistic values that 
agree to within 1%. Since each file contains the 13 000 
events with the largest values of the ^-statistic, numer- 
ical fluctuations in determining the value of T can lead 
to slightly different lists being returned on different plat- 
forms. Hence the validator tolerates unpaired candidate 
events if they lie within 1% of the smallest value on the 
list. 

A workunit is validated once it has three results that 
agree with one another to within these tolerances. If 
the three results do not pass this validation step, the 
Einstein@Home server generates more instances of this 
workunit until three valid results have been obtained. 
For the search described in this paper, a "post-mortem" 
analysis of the computation shows that the probability of 
a successful but invalid result is small (0.36%), and the 
errors which make a successful result invalid are typically 
unique and irreproducible. Hence we estimate that it is 
highly improbable that even a single incorrect result has 
been marked as "valid" by the automatic validator. 



B. Workunit design 

The different workunit s cover (search) different parts 
of parameter space. A key design goal of these workunit s 
is that they should have roughly equal computational run 
times, of the order of ~ 8 hours, and that the computa- 
tional effort to carry out the entire search should last 
about 0.5-1 years. Another key design goal is that each 
workunit uses only a small re-usable subset of the total 
data set. These allow Einstein@Home volunteers to do 
useful computations on the one-day timescale, and mini- 
mizes the download burden on their internet connections 
and on the Einstein@Home data servers. 

Each workunit uses only one segment of data over a 
narrow frequency-range, but covers the whole sky and the 
full range of frequency- derivatives / at that frequency. 
Therefore, the entire search is divided into computational 
units over different data segments and frequency-bands. 
In the following it will be useful to label the workunits 
by three indices (j, /c, £), where j = 1, • • • , 17 denotes the 
data segment, /c = 1, • • • , 2900 labels the 0.5 Hz band cov- 
ered by the input data file, and ^ = 1, • • • , M(j, /c) enu- 
merates individual workunits associated with data seg- 
ment 2 and frequency band k. Note that each workunit 
uses a frequency band that is smaller than the 0.5 Hz 
covered by the input data files, i.e. M(j, k) > 1. 



1. Search grid in parameter space 

The parameter space is gridded in such a way that 
no point has a "squared-distance" from its nearest grid 
point that exceeds a certain "maximal mismatch". The 
distance is defined by a metric on parameter space, first 
introduced in [2TJ [22] . The squared distance is the frac- 
tional loss of squared signal-to-noise ratio (SNR ) due 



to waveform mismatch between the putative signal and 
the template. The search grid was constructed based 
on the projected metric on the subspace orthogonal to 
the frequency direction df with / = 0.. For any given 
workunit, the parameter-space grid is a Cartesian prod- 
uct of uniformly-spaced steps df in frequency, uniformly- 
spaced steps df in frequency derivative, and a two- 
dimensional sky-grid, which has non-uniform spacings 
determined by the projected metric. For frequencies in 
the range [50, 290) Hz, the maximal mismatch was chosen 
as 771 = 0.2 (corresponding to a maximal loss in SNR of 
20%), while in the range [300, 1500) Hz, the maximal mis- 
match was m = 0.5. Due to a bug in the script generating 
the sky-grids, the range [290, 300) Hz, was covered by fre- 
quency and spin-down steps corresponding to m = 0.2, 
whereas the sky-grids were constructed for m = 0.5. The 
distribution of actual mismatches in this frequency range 
will therefore be somewhat in between those of the low- 
frequency and high-frequency workunits. 

It can be shown [45 that these relatively large mis- 
matches give near-optimal sensitivity for a coherent 
search at fixed CPU power. Choosing finer grid spac- 
ings (i.e. a smaller mismatch) would require searching 
more grid-points, thus reducing the maximal possible co- 
herent integration time. A coarser search grid would al- 
low longer integrations but at a larger average loss in 
SNR. Because of these two competing tendencies, the 
sensitivity as a function of mismatch m has a maximum 
in the range m ~ 0.25-0.7, depending on the choice of 
false-dismissal rate from the grid mismatch. Full details 
of the parameter-space grid and workunit construction 
are given in [45 ; a short summary follows. 



2. Search grid in frequency and frequency- derivative 

The step-size in frequency was determined using the 
metric-based expression 



dfj = 



■kT, 



(1) 



span, J 



SO the frequency-spacing depends on T^panj of the 
data segment j. For the low- frequency workunits 
(/ < 300 Hz), this results in frequency steps in the range 
dfj G [3.43, 4.06] /iHz, while for high-frequency workunits 
dfj e [5.42, 6.41] /iHz. 

The range of frequency-derivatives / searched is de- 
fined in terms of the "spin-down age" r = — ///, namely 
r > 1000 years for low-frequency and r > 10 000 years for 
high-frequency workunits. Younger neutron stars than 
the limited range of the search probably would have left 
a highly visible (Sedov phase) supernova remnant or a 
pulsar wind nebula, and thus our search for unknown 
neutron stars targeted older objects, which also resulted 
in less computational cost. The search also covers a 
small "spin-up" range, so the actual ranges searched are 
/ G [— //t, 0.1//r]. In / the grid points are spaced ac- 
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each workunit is chosen to be approximately constant for 
ah workunits. However, in practice, this number can vary 
by up to a factor of 2 due to discretization effects. The 
number of points in the sky-grids grows with frequency as 
/^ and the number of points in the spin-down grid grows 
hnearly with /. Thus, to keep the number of templates 
(and therefore the CPU time) approximately constant, 
the frequency range covered by each workunit decreases 
as f~^. Hence for fixed j, M(j, k) is roughly proportional 
to k^. 



C. The output of a workunit 



FIG. 3: Four different sky-grids in Hammer- Aitoff |46] projec- 
tion. The top row is for frequency / = 60 Hz and the bottom 
row is for / = 310 Hz. The left column shows data segment 
j = 1 (from HI) with a spanned time of Tspan,i = 33.8 h, 
while the right column shows data segment j = 15 (from LI) 
with a spanned time of Tspan,i5 = 39.8 h. The grid points 
are spaced more closely for a longer spanned time, and for a 
higher frequency. 



cording to 



df: 



^ 7rT2 

span, J 



(2) 



resulting in resolutions dfj G [1.84, 2.59] x 10"^°Hz/s 
for low- frequency workunits, and dfj G [2.92, 4.09] x 
10~^^ Hz/s for high-frequency workunits, depending on 
the duration T^panj of different segments j. 



3. Search grid in the sky parameters 



The resolution of the search grid in the sky-directions 
depends both on the start-time tj and duration Tspan,j of 
the segment, as well as on the frequency. The number of 
grid points on the sky scales as ex /^ , and approximately 

as (X ^span,j fo^ ^h^ range of Tspan,j ~ 30 40 h used 

in this search. Contrary to the simple uniform spacings 
in / and /, the sky- grids are computed beforehand and 
shipped together with the workunits. In order to sim- 
plify the construction of workunits and limit the amount 
of different input-files to be sent, the sky-grids are fixed 
over a frequency range of 10 Hz, but differ for each data 
segment j. The sky-grids are computed at the higher end 
of each 10 Hz band, so they are slightly "over- covering" 
the sky instead of being too coarse. The search covers 
a frequency range of 1450 Hz, and so there are 145 dif- 
ferent sky- grids for each segment. To illustrate this, four 
of these sky-grids are shown in Figure |3] corresponding 
to two different data segments at two distinct frequency 
bands. 

To ensure that each workunit takes a similar amount 
of CPU time, the total number of template grid points of 



The result from completing one workunit on an Ein- 
stein@Home host computer is a ZIP-compressed ASCH 
text file containing the 13 000 candidate events with the 
largest values of the ^-statistic found over the parameter- 
space grid points analyzed by that workunit. Each line 
of the output file contains five columns: frequency (Hz), 
right ascension angle (radians), declination angle (ra- 
dians), spin-down-rate (Hz/s) and 2jr (dimensionless). 
The frequency is the frequency at the Solar System 
Barycenter (SSB) at the instant of the first data point 
in the corresponding data segment. 

The number 13 000 was decided in advance, when the 
workunits were first launched on the Einstein@Home 
project, which was about one year before the post- 
processing pipeline was developed. The network band- 
width required to retain more than 13 000 candidates 
per workunit, and the storage space required to pre- 
serve them, would have exceeded the capacity of the 
Einstein@Home server and its internet connection. For 
frequency-band and data segment combinations with 
small numbers of workunits, for example the j = 1 data 
set from 301.0 Hz to 301.5 Hz, almost ah of the 13 000 
candidate events are later used in the post-processing 
pipeline. However (as can be seen later in Figure B for 
most frequency-bands the post-processing pipeline only 
needed and used a fraction of the events that were re- 
turned. 

Returning the "loudest" 13 000 candidate events effec- 
tively corresponds to a fioating threshold on the value of 
the ^-statistic. This avoids large lists of candidate events 
being produced in regions of parameter space containing 
non- Gaussian noise, such as instrumental artifacts that 
were not removed a priori from the input data. 



IV. POST-PROCESSING 

As shown previously in Figure [2J after result files are 
returned to the Einstein@Home servers by project par- 
ticipants, further post-processing is carried out on those 
servers and on dedicated computing clusters. The goal 
of this post-processing pipeline is to identify consistent 
candidate events that appear in many of the 17 different 
data segments. 



In this paper, a consistent (coincident) set of "candi- 
date events" is called a "candidate". Candidate events 
from different data segments are considered coincident 
if they cluster closely together in the four-dimensional 
parameter space. A clustering method using a grid of 
"coincidence cells" will reliably detect strong CW sig- 
nals, which would produce candidate events with closely- 
matched parameters. 

The post-processing pipeline operates in 0.5 Hz- wide 
frequency-bands, and can be summarized in three steps. 
In step one^ the coincident false alarm probability is fixed. 
In step two^ the frequency values of candidate events are 
shifted to the same fiducial time. In step three, a grid 
of cells is constructed in the four-dimensional parameter 
space, and each candidate event is assigned to a particu- 
lar cell. In the following the details involved in each step 
are described. 
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A. Preparation and selection of candidate events 

In the first step the individual result files are prepared 
for the later analysis by uncompressing them and keeping 
only a subset of the candidate events: from the (j, /c, ^)'th 
workunit only the S{j,k,£) candidate events with the 
largest values of 2jr are retained. 

The number of these candidate events is chosen a priori 
to obtain a pre-determined fixed false alarm probability. 
The false alarms should be approximately uniformly dis- 
tributed among the workunits, since each workunit exam- 
ines a similar number of independent grid points in pa- 
rameter space. The number of candidate events is chosen 
so that in a 0.5 Hz- wide frequency-band the probability 
that one or more coincidence cells after doing the clus- 
tering (in step three) has Cmax = 7 or more coincidences 
is Pf = 0.001. Thus, in the analysis of 2900 such fre- 
quency bands, in random noise one would expect to find 
only about three candidates with sev en or m ore coinci- 
dences. (As explained later in Section IVF 1 , this overall 



probability for the entire search is somewhat increased 
because the coincidence cell grids are also shifted by half 
their width in 16 possible combinations). 

In terms of the notation introduced in the previous 
section, the number of candidate events kept from the 
(j, k, iy th workunit takes the form 



£{j,k,e) 



M{j,k) 



(3) 



where Sseg{k) is shown in Figure |4J Because the indi- 
vidual workunits are constructed to use approximately 
the same amount of CPU time, each workunit examines 
approximately the same number of templates in param- 
eter space, so the same number of candidate events are 
retained from all workunits which have the same input 
data file and data segment. This implies that the number 
of candidate events that are kept per data segment j and 



FIG. 4: The circles show the number of candidate events 
Sseg{k) retained per data segment and per 0.5 Hz frequency 
band in the post-processing in each 10 Hz band. The 
dashed curve represents the number of candidate events which 
are returned from volunteering hosts participating in Ein- 
stein@Home. The strange location of the point at 290 Hz 
is explained in Sec. |IIIB 1[ 



per frequency band is independent of the data segment j: 



M{j,k) 



,{k) 



(4) 



Since the sky-grids are fixed in 10 Hz intervals, £seg{k) 
takes the same value for all values of k in the range of 
20p + 1, • • • , 20(p + 1) where p labels the sky-grids by an 
integer in the range p G 0, • • • , 144. 

It is illustrative to look at a specific case. For exam- 
ple consider the 0.5 Hz band covering [301.0, 301.5) Hz, 
this band is labeled hy k = 503. As is shown in Fig- 
ure [4j in this band the post-processing pipeline retains 
^seg(^ = 503) = 24 960 candidate events from each of the 
17 different 30- hour data segments. The 30-hour data 
segment from HI with the shortest time span (j = 1) has 
approximately 4.3 x 10^ templates divided among just 
M{j = l,k = 503) = 2 workunits, so 12 480 candidate 
events are retained from each of these workunits. The 30- 
hour data segment from LI with the longest time span 
{j = 15) has approximately 1.7 x 10^ templates divided 
among M{j = 15,/c = 503) = 7 workunits, so 3 565 can- 
didate events are retained from each of these workunits. 
In the later stage of the post-processing, this ensures that 
each of the different data segments contributes equally to 
the probability of generating false alarms in the coinci- 
dence step. 



B. Number of cells in the post-processing 
coincidence grid 

It is important to calculate the number of coincidence 
cells in the coincidence grid. Together with the desired 
false alarm probability, this determines the number of 
candidate events to retain in the post-processing pipeline. 

The number of coincidence cells Ncei\{k) contained in 
each 0.5 Hz frequency band k (while doing the clustering 
in step three) is determined by the sizes of the cells. This 
is given by 
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where A/ denotes the coincidence cell width in frequency, 
A/ denotes the width in spin-down, and Aa{S) and 
AS{S) denote the coincidence cell widths in right ascen- 
sion and declination (both of which vary with declination 
5). The choice of the coincidence cell sizes will be ex- 
plained in detail later when step three will be described. 



In random noise, the candidate events produced by each 
different data segment are independent, so that the coin- 
cidence cells that are "marked" by one or more candidate 
events are also described by a (different) binomial distri- 
bution. Without loss of generality, again consider the 
first coincidence cell. The probability that it contains 
candidate events from n distinct data segments is then 
given by 
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is the binomial coefficient. Thus the 



where (^) = ^^_^y^ 

probability per coincidence cell of finding Cmax or more 

coincidences is given by 
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For the desired Pf = 0.1% = 10"^ per 0.5 Hz band k, 
this equation is solved numerically to find S^egik). The 
results for Sseg{k) are shown in Figured 



False alarm rate and the number of candidate 
events retained 



D. Choice of false alarm probability and detection 
threshold 



The number of candidate events that must be retained 
is determined by the number of cells in the coincidence 
grid A^ceii(^) and by the desired probability of false alarm 
Pp for false coincidence of candidate events from Cmax 
or more data segments in each 0.5 Hz band. To re- 
late these quantities, consider the case of random in- 
strumental noise, in which the candidate events are dis- 
tributed uniformly about the coincidence grid. Concen- 
trate on a single 0.5 Hz band /c, and consider the first 
of the A^seg = 17 data segments. A total of Eseg{k) 
candidate events must be distributed uniformly among 
Nceii{k) coincidence cells. Each candidate event falls in 
a random coincidence cell, independent of the locations 
of the previous candidate events. The probability that 
the first candidate event falls in the first coincidence cell 
is 1/Nceii{k)^ and hence the probability that the first co- 
incidence cell remains empty is 1 — 1/A^ceii(^)- If the re- 
maining Eseg{k) — 1 candidate events fall independently 
and at random into the coincidence cells, then this gen- 
erates a binomial distribution, and the probability that 
the first coincidence cell contains no candidate events is 
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Since the first coincidence cell is equivalent to any other, 
the probability that the candidate events from the first 
data segment populate any given coincidence cell with 
one or more candidate events is thus given by 



e{k) = l-pk{0) = l 
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The goal of this work is to make a confident detec- 
tion^ not to set upper limits with the broadest possible 
coverage band. This is reflected in the choice of the ex- 
pected false alarm probability and the choice of a detection 
threshold. 

The detection threshold of 12 events was chosen be- 
cause, as described in Section |VII[ the hardware injec- 
tions are only "turned on" for 12 of the 17 data seg- 
ments. The detection threshold ensures that these simu- 
lated signals are properly detected by the post-processing 
pipeline. 

The choice of false alarm probability (Pp = 0.1% = 
10~^ per 0.5 Hz band to have coincidences in Cmax = 7 
or more data segments) is a pragmatic choice, which leads 
to an extremely small false alarm rate at the detection 
threshold. For actual data, the probability of finding 
7 or more coincidences in a given 0.5 Hz band can be 
somewhat larger than the target value of 0.1% because 
the candidate events are not uniformly distributed over 
the grid of coin cidence cells and because (as described in 
Section IVFl) 16 sets of coincidence cells are used for 



each 0.5 Hz band. 

In random noise, the probability of reaching the detec- 
tion threshold of 12 coincidences depends on the number 
of cells in the coincidence grid, which is a function of 
frequency. Some representative numbers are given in Ta- 
ble HI they vary from about 10~^^ to 10~^^ depending 



upon the 0.5 Hz band. The false alarm probabilities de- 
crease very rapidly with increasing coincidence number. 
For example the probability of finding 14 or more co- 



incidences in random noise varies from about 10 
10-2^ 
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to 
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/fc[Hz] k iVceii(A;) Pf{C > 7) Pf{C > 12) Pf{C > 14) 

1.5 X 10"^^ 3.0 X 10"^^ 
8.7 X 10"^^ 5.7 X 10"^° 
6.7 X 10"^^ 9.9 X 10"^^ 
2.2 X 10"^^ 8.4 X 10"^^ 



50.0 1 734 500 10" 

290.0 481 35198 800 10" 

301.0 503 2161284 10" 

1499.5 2900 233 345 981 10" 



299.01 - 1 - 299.51 



TABLE III: False alarm probabilities Pp in four different 
0.5 Hz frequency bands labeled by the integer k. The fre- 
quency at the lower boundary of the 0.5 Hz band k is denoted 
by /fc. The number of coincidence cells in the /c'th half- Hz fre- 
quency band is denoted by Nceiiik). The probability of find- 
ing 7 or more coincidences {C > 7) in randomly-distributed 
noise is fixed to be 0.1%. The probability of finding 12 or 
more coincidences (the detection threshold, C > 12) in ran- 
dom noise varies over two orders of magnitude, from about 
10~^^ to 10~^^. The probability of finding 14 or more coinci- 
dences {C > 14) in random noise varies from about 10~^^ to 

io-^\ 



Once might ask why we chose to specify a uniform 
false alarm probability, across all frequencies, of 0.1% for 
Cmax = 7, rather than directly specify a much lower false 
alarm probability at the detection threshold C = 12. This 
was because we wanted the most significant coincident 
events due to noise alone to have C values a few less than 
our detection threshold, and we wanted these candidates 
to be uniformly distributed over frequency bands. Any 
detected signals could then be compared against fairly 
uniform fields of noise candidates in adjacent frequency 
bands. If a uniform false alarm probability had been 
specified at the C = 12 level, then the expected noise 
candidates with C ^ 7 would not have been uniformly 
distributed over frequency, due to the differing numbers 
of coincidence cells in each frequency band. 

The choice of detection threshold and false alarm prob- 
ability sacrifice a small amount of sensitivity compared 
with a higher values, but ensures that high numbers of 
coincidences are extremely improbable in random noise. 
A strong signal (say a factor of 1.5 above the upper curve 
in Figure [9| would be expected to produce 15 or more co- 
incidences in this detection pipeline. With the thresholds 
that we have adopted, this would stand out very strongly: 
the probability of having even one such an event appear 
in coincidence in random noise is about 10~^^ per 0.5 Hz 
band. 



E. Shifting candidate event frequencies to a fixed 
fiducial time 



In the second step of the post-processing, the frequency 
value of each retained candidate event is shifted to the 
same fiducial time: the GPS start time of the earliest 
{j = 4) data segment, tfiduciai = t^ = 793 555 944 s. This 
shifting is needed because a CW source with non-zero 
spin-down would produce candidate events with different 
apparent frequency values in each data segment. This 
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FIG. 5: Additional "wings" at the boundaries of each 0.5 Hz 
frequency band must be included in the coincidence analy- 
sis stage of the post-processing. This is because spin-down 
can carry a source below this half-Hz band, and spin-up can 
carry it above the band. To illustrate this, the frequency 
band k = 498 (covering [299, 299.5) Hz) is (partly) shown by 
the dark-gray shaded area. The dashed sloped lines show the 
boundaries of the small additional regions (light gray) in fre- 
quency space whose candidate events must also be considered 
in the post-processing. Because the allowed spin-up range is 
ten times smaller than the allowed spin-down range, the up- 
per boundary has a slope ten times smaller than the lower 
boundary. 



step would shift these candidate events back to the same 
frequency value: 



/(^fiducial) = f{tj) + [^fiducial 



tj] f , 



(10) 



where / and f{tj) are the spin-down-rate and frequency 
of a candidate event reported by the search code in the 
result file, and tj is the timestamp of the first datum 
in the data segment. At the end of the second step, all 
candidate events for the 0.5 Hz band are merged into one 
file. 

These candidate events are collected from a frequency 
interval that is slightly wider than 0.5 Hz. To see why 
this is necessary, consider a potential source whose fre- 
quency in the first data segment (j = 4) is at the lower 
(or upper) boundary of the 0.5 Hz interval. If the source 
has the minimum (or maximum) allowed value of /, then 
in the later data segments it moves into, or is recorded 
in, the previous (or next) 0.5 Hz band. This effect is 
most apparent for the last j = 11 data segment, as illus- 
trated in Figure [5| So in collecting the candidate events 
for analysis of a given 0.5 Hz band, the frequency range 
is enlarged slightly for events coming from later and later 
data segments, as shown in Figure [5] 
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F. Search for coincident candidate events 

The third step and final stage of the post-processing 
is to search for parameter-space coincidence among the 
candidate events. If a CW source is present that is strong 
enough to be confidently detected, then it would produce 
large JT-statistic values (i.e. candidate events) in many 
or all of the 17 data segments. In addition, the values of 
the frequency at the fiducial time /(tfiduciai), sky position 
(given by right ascension a and declination 5), and spin- 
down / for these candidate events would agree among 
all data segments (within some coincidence "window" or 
"cell"). 



1. Coincidence search algorithm 

To find coincidences, a grid of cells is constructed in 
four-dimensional parameter space, as described previ- 
ously. This analysis uses rectangular cells in the coor- 
dinates (/, /, a cos ^, 6). The dimensions of the cells are 
adapted to the parameter space search. Each candidate 
event is assigned to a particular cell. In cases where two 
or more candidate events from the same data segment j 
fall into the same cell, only the candidate event having 
the largest value of 2J^ is retained in the cell. Then the 
number of candidate events per cell coming from distinct 
data segments is counted, to identify cells with more co- 
incidences than would be expected by random chance. 

The search for coincident cells containing large num- 
bers of candidate events is done with an efficient code 
that uses linked-list data structures, A^logA^ sort algo- 
rithms, and log A/" bisection search algorithms. To ensure 
that candidate events located on opposite sides of a cell 
border are not missed, the entire cell coincidence grid is 
shifted by half a cell width in all possible 2^ = 16 combi- 
nations of the four parameter-space dimensions. Hence 
16 different coincidence cell grids are used in the analysis. 

The cells in the coincidence grid are constructed to 
be as small as possible to reduce the probability of co- 
incidences due to false alarms. However, since each of 
the 17 different data segments uses a different parame- 
ter space grid, the coincidence cells must be chosen to 
be large enough that the candidate events from a CW 
source (which would appear at slightly different points in 
parameter space in each of the 17 data segments) would 
still lie in the same coincidence cell. 



2. Frequency and spin-down coincidence windows 

In frequency, the spacing of the parameter-space grid 
is largest for the data segment with the smallest value 
of Tspan,j, which is the first data segment j = 1. At 
first, this would appear to be the correct size A/ for 
the coincidence cell in the frequency direction. However 
since the frequency of a candidate event must be shifted 
to a fixed fiducial time according to its spin-down value. 



and since that spin-down value can not be more accurate 
than the / spacing, the size of the coincidence cell must 
be somewhat larger to accommodate the effects of this 
discretization error in /. The coincidence window in the 
frequency direction is thus determined by 



A/ = max (dfj + At dfj 



(11) 



where the maximization over j selects the data segment 
with the smallest T^panj (which is j = 1) and 

At = I max tj - min tj \ = tn - t^ =1 997 256 s (12) 
j j 

is the total time span between the latest and earliest data 
segments. For safety, e.g. against noise fluctuations that 
could shift a candidate peak. A/ has been increased by 
a further 40% below 300 Hz, so that the width of the 
coincidence cell in frequency is A/ = 0.77 x 10~^ Hz, 
and by 30% above 300 Hz, so that A/ = 1.08 x 10"^ Hz. 
For the spin-down parameter, the size of the coinci- 
dence cell is given by the largest dfj spacing in the pa- 
rameter space grid, which is also determined by the small- 
est value of Tspan,j- For safety this is also increased by 
40% below 300 Hz giving A/ = 3.7 x 10"^° Hz s"\ and 
by 30% above 300 Hz giving A/ = 5.18 x 10"^° Hz s"^ 



3. Coincidence windows in apparent sky position 

Determining the optimal size for the coincidence cells 
in the sky coordinate directions is more difficult. Each of 
the 17 different data segments uses a different sky- grid, as 
illustrated in Figure [3J Ideally the size of the coincidence 
cells in these sky directions must be just large enough to 
enclose one parameter space grid point from each of the 
17 different sky- grids. A practical solution to determine 
the coincidence cells, which is close to optimal, makes use 
of an observation concerning the parameter-space metric 
that first appears in [25 . 

To understand the properties of the parameter-space 
metric, first consider the relative orders-of-magnitude of 
the different frequency modulation effects. The frac- 
tional Doppler shift due to the Earth's annual orbital 
motion about the Sun has magnitude |vorbitai|/c = 10~^ 
and the fractional Doppler shift due to the detector's 
daily motion about the Earth rotation axis has magni- 
tude |vrotationai|/c = 10"^. For the Tspan ^ 40 h period 
of a single coherent integration, one can approximate the 
motion of the Earth's center of mass as motion with con- 
stant acceleration (along a parabolic trajectory) rather 
than as circular motion. The neglected term in the frac- 
tional Doppler shift has magnitude |vorbitai|^span/2^ ^ 
|vorbitai||^pT2p^n/2c ~ 4 X 10-^ whcrc |u;| = 27r/year 
is the magnitude of the Earth's orbital angular velocity 
about the sun. This term is a factor of 25 smaller than 
|vrotationai|/c and hcucc cau be neglected. With this ap- 
proximation, the orbital motion of the Earth is simply 
responsible for an apparent shift in the frequency / and 
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FIG. 6: Example sky- grid and its projection onto the equa- 
torial plane. This sky- grid corresponds to the data segment 
j = 7 used in the frequency range from 60 Hz to 70 Hz. The 
top plot shows a Hammer-Aitoff projection of the sky-grid. 
The middle plots show the projection of the sky-grid points 
in the northern hemisphere (left column) and in the southern 
hemisphere (right column) onto the equatorial plane. The 
bottom plots show histograms of cos((5) and the dashed line 
represents a linear fit to the distribution showing its unifor- 
mity. 



spin-down rate / of a source: the effects of the Earth's 
center of mass motion are degenerate with a shift in fre- 
quency and spin-down. So the Earth's orbital motion 
causes a signal only to shift to a different template in / 
and /; the Earth's rotation has a period of one sidereal 
day and can not be modeled by a shift in / or /. Note 
that terms are neglected only in determining where to 
place search grid points in parameter space (because the 
neglected terms have an insignificant effect on where the 
grid points are placed). The actual filtering of the data 
uses "exact" barycentering routines (total timing errors 
below 3/is). 



The search grid in parameter space is a Cartesian prod- 
uct of a frequency grid, a spin-down grid, and a two- 
dimensional sky grid. Since the search maximizes the de- 
tection statistic over frequency and spin-down, the met- 
ric used to place grid points on the sky [45 may be ob- 
tained by minimizing the four-dimensional metric over 
frequency and spin-down and projecting it into the sky 
directions. As shown in the previous paragraph, this two- 
dimensional projected sky metric is well- approximated by 
assuming that the Earth is spinning about its axis but 
has its center of mass at rest. If the coherent integration 
period is an integer number of days, then by symmetry 
the two-dimensional metric on the sky is invariant under 
rotation about Earth's axis (pa. is a Killing vector). This 
is still an approximate symmetry for the search described 
here, since the coherent integration period and Tgpan are 
longer than the rotation period (one day). 

One can easily find coordinates in which this approxi- 
mate sky metric (the four-dimensional metric, minimized 
over frequency and spin-down and projected onto the sky 
directions) is proportional to diag(l, 1). These new sky- 
coordinates are obtained by perpendicular projection of a 
point on the two-sphere (celestial sphere) vertically down 
into the unit radius disk that lies in the equatorial plane. 
If n denotes a unit vector pointing from the SSB to the 
source the new coordinates are the components of n in 
the equatorial plane: Ux — cos J cos a, Uy = cosS sin a. 
Points which are equally spaced in these coordinates cor- 
respond to equal spacing in Doppler shift, since source 
Doppler shift due to the Earth's rotation is just pro- 
portional to the component of the source direction vec- 
tor in the equatorial plane. It then follows from rota- 
tional invariance that (with these approximations) the 
projected sky metric in these coordinates is proportional 
to diag(l,l) [43. The effect may be immediately seen 
in Figure |6j the grid of "equally-spaced" points forms a 
(roughly) uniform square grid on the unit radius disk in 
the equatorial plane. Computing the Jacobian of the co- 
ordinate transformation shows that in the original coor- 
dinates (a, S) the coordinate-space density of grid points 
should be proportional to | cosS smS\ = | sin(2(5)|. 

This simple behavior of the projected sky metric guides 
the construction of the coincidence-windows in the sky 
directions. Define polar coordinates (r, a) on the unit 
radius disk in the equatorial plane by r = cos S. The co- 
ordinate boundaries of uniformly distributed coincidence 
cells containing a single parameter-grid point are then 
given hy r da = dr = const. When written in terms of 
the original coordinates this becomes 



cos(^) da = I sin((5)| dS = const. 



(13) 



This is not directly useful, because it is singular as 5 ^ 0, 
but suggests a coincidence window size which varies with 
declination according to the model 



Aa{5) 
A5{5) 



Aa{0)/cos{5) 

i A(5(0) 

1 Aa(0)/|sin( 



(14) 

if |<^| < S„ 
KAa{0))\ if \5\>S^. 
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To ensure continuity at ^ = ^c, the transition point Sc is 
defined by the condition AQf(0)/| sin(|Jc| — i^Aa{0))\ = 
A^(0). Ai: is a tuning parameter of order unity, described 
below. An example of this coincidence window model is 
shown in Figure [7[ 

For each of the 145 different 10 Hz bands, the win- 
dow size is determined by the three constants Aa(0), 
A(5(0) and hi. For each sky-grid p these values are di- 
rectly determined from the sky-grids used in the search 
as follows. For each 10 Hz frequency band the maximum 
distances between adjacent declination points to either 
side are calculated for each of the 17 sky- grids as a func- 
tion of declination 6. In this way, 17 different overlaying 
curves Aj{S) (one per data segment) are obtained. These 
are indicated by the circles in Figure [7| for a representa- 
tive 510 Hz — 520 Hz frequency band as illustration. Then 
the parameter A(5(0) is obtained by considering the max- 
imum separation to either side between all neighboring 
declination grid points and between the 17 different sky- 
grids, increased by a 10% safety factor as 



AS{0) = 1.1 max{A^((5)} . 



(15) 



The largest separations near the poles (1.4 < \S\ < 7r/2) 
are then found and increased by a safety factor of 20% 
to determine the parameter Aa{0) via 



Aa(0) 



1.2 max{A^( 

j,5 



> 1-4)} . 



(16) 



Finally, the parameter k, was chosen by visually examin- 
ing diagrams similar to Figure [7| for all 145 of the 10 Hz 
bands. A k. value of 1.5 was found to be sufficient in 
most cases, while some bands required somewhat higher 
or lower values. For each triple of sky-coincidence param- 
eters, tests were then performed to check that each sky- 
cell contained at least one sky-point from each data seg- 
ment. In Figure [T] the complete declination coincidence 
window model given by Equation (14) is represented by 
the solid black curve. 

The three parameters for all sky-grids as a function of 
frequency are shown in Figure [SJ As stated above, the 
sky-grids are constant for 10 Hz- wide steps in frequency, 
and so these parameters vary with the same step-size. 



G. Output of the coincidence search and 
significance of a candidate 
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FIG. 7: The sky coincidence- window model for the frequency 
band from 510 — 520 Hz. The horizontal axis shows the dec- 
lination S in radians. On the vertical axis, the circles labeled 
Aj (6) correspond to the the maximum distance in radians to 
neighboring ^-points on either side. The solid curve shows the 
declination coincidence- window model A^(^) with the param- 
eters A^(0) = 0.2489, Aa(0) = 0.0433, and k = 1.5 used in 
this frequency band. It lies just above the largest declination 
separations shown. The stars denote the borders of the dec- 
lination coincidence cell-grid and the diamonds represent the 
borders of the shifted declination coincidence cell- grid. 



0.8 



"O 

x:o.6 

to 

<1 
§0.4 

< 



0.2 



' 


' ' 


' 




X Aa(0) 




f 


+ A5(0) 










+ 




+ 




1 






i + 


1. 




V ^HK; 


ii^-i- 


w: 




m 



2.5 



Q. 
Q. 
03 



0.5 



50 500 1000 1500 
Frequency [Hz] 



50 500 1000 1500 
Frequency [Hz] 



The output of the post-processing pipeline is a list 
of the most coincident candidates. In each frequency 
band of coincidence- window width A/, the coincidence 
cell containing the largest number of candidate events is 
found. Thus for each frequency band the pipeline finds 
the most coincident candidate maximized over the sky 
and over the spin-down parameter-range. The pipeline 
outputs the average frequency of the coincidence cell, 
the average sky position and spin-down of the candidate 
events, the number of candidate events in the coincidence 
cell, and the significance of the candidate. 



FIG. 8: The parameters Aa(0), A(5(0) and k of the sky 
coincidence- window model as a function of the 10 Hz fre- 
quency band. 



The "significance" of a candidate was first introduced 
in [27 . A candidate, consisting of the candidate events 
1, . . . , Q, has significance 



S{1,...,Q) 



Q 

5=1 



ln(l + ^J) 



(17) 
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where Q < 17 is the number of candidate events in the 
same coincidence cell. To understand the meaning of 
the significance, consider the case of pure Gaussian noise 
with no signal present. In this case the values of 2J^ have 
a central x^ distribution with four degrees of freedom. 
The corresponding probability density function po of 2jr 
is given by 



Po(2^) = ^ e- 



(18) 



The false alarm probability Pq that 2J^ exceeds a certain 
threshold 2J^q when there is no signal present has the 
form 



Po(2^o) = (l + ^o)e' 



-^0 



(19) 



The joint false alarm probability of candidate events 
1, . . . , Q can be written as flgLi ^o(2J^g). Therefore, in 
this analysis candidates are ranked according to 



n^o(2^,) 



(20) 



where S = Ylq=i ~l^^o(2J^g) is exactly the significance 



defined in Equation (17). Thus ranking candidates by S 



is equivalent to ranking them by false alarm probability: 
candidates with large positive significance would not be 
expected to occur in Gaussian random noise. As will be 



described later in Section VIII the significance is used 
to rank equally coincident candidates within the same 
narrow frequency-band. In such cases the candidate with 
the largest significance is considered. 

The post-processing pipeline has been validated by in- 
ternal testing, and also using simulated CW signals cre- 
ated via so-called "software injections" [48 . In addition. 
Section |VII| presents realistic end-to-end testing of the 
analysis pipeline using "hardware injections" , where sim- 
ulated isolated-pulsar signals are physically added into 
the interferometer control systems to produce instru- 
mental signals that are indistinguishable from those that 
would be produced by physical sources of gravitational 
waves. 



V. ESTIMATED SENSITIVITY 



is of order 10"^"^, making it extremely unlikely for can- 
didate events from random noise to show up consistently 
in 12 or more segments of data. This is therefore an esti- 
mate of the signal strength required for high-confidence 
detection. The pipeline developed for this purpose oper- 
ates in 0.5 Hz frequency bands and consists of testing a 
large number of distinct simulated sources (trials) to see 
if they are detectable. A "trial" denotes a single simu- 
lated source which is probed for detection. 



A. Iteration method 

For every trial, source parameters are randomly cho- 
sen independent of the previous trial, except for the in- 
trinsic amplitude hp. For the first trial ho is set to a 
starting value 30^s/Sh/^0 hours. The rule for varying ho 
depends upon the last A/'iast trials, where N^^^"" = 100, 
^\fsf = 20, and N^^^f = 100. In the past Mast trials, if 
more than 10%, 50%, or 90% of simulated sources have 
been detected then ho is decreased by 0.25 ho/ntria\ for 
the following trial, where ntriai is an integer in the range 
< ntriai < 1000 that is incremented with each addi- 
tional trial. On the other hand, if less than 10%, 50%, 
90% of simulated sources have been detected then ho is 
increasedhy 0.25 /^-o/^triai for the next trial. This process 
is followed until ho has converged to a stationary range 
after 1000 trials. Then the median of ho is found using 
the /iQ-values starting from that trial, where the desired 
detection level has been reached the first time during the 
A/'iast trials. The following describes the pipeline for a 
single trial. 



B. Population of simulated sources 

For each trial, a random source frequency is chosen 
from a uniform distribution in a given 0.5 Hz frequency 
band and a frequency-derivative is drawn from a uniform 
distribution in the range covered by the Einstein@Home 
search. A sky position is chosen from a uniform distri- 
bution on the celestial sphere, and random values are 
chosen for the pulsar "nuisance parameters". These are 
the inclination parameter cos(^), initial phase ^o? ^ind 
polarization angle i/j as defined in [19] , and are all drawn 
from the appropriate uniform distributions. 



The sensitivity of the search is estimated using Monte- 
Carlo methods for a population of simulated sources. The 
goal is to find the strain amplitude ho at which 10%, 
50%, or 90% of sources uniformly populated over the 
sky and in their "nuisance parameters" (described below) 
would be detected. As previously discussed, the false 
alarm probability (of getting 7 or more coincidences in a 
0.5 Hz frequency-band) is of order 10~^. In this analysis, 
"detectable" means "produces coincident events in 12 or 
more distinct data segments" . The false alarm probabil- 
ity for obtaining 12 or more coincidences in a 0.5 Hz band 



C. Determination of 2J^ values for a single 
simulated source 



The noise fioors of the different SFTs are estimated at 
the source's frequency intervals using a running median 
with a window of ±25 frequency bins. Figure [l] showed 
the average of these for the data segments used from each 
instrument. 

Then for each set of parameters the detection statis- 
tic 2J^ is estimated using a semi-analytic method. From 
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the estimated noise floor at the simulated source's fre- 
quency and given the other source parameters, the expec- 
tation value of the ^-statistic is calculated analytically 
as given in [W . A random number is then drawn from a 
non-central x^ distribution with four degrees of freedom 
and with the corresponding mean value. 

These random numbers, drawn from the appropriate 
distribution of 2J^ values, would be sufficient to deter- 
mine the sensitivity of the search, if the template grid 
in parameter space were very closely spaced, so that the 
template bank always contained at least one waveform 
that was a very close match to the putative signal. How- 
ever the grid in parameter space used in this search is 
quite "coarse" , corresponding to a mismatch of 20% be- 
low 290 Hz and 50% above 300 Hz, so that the 2J=' value 
that would be returned by the search might be signifi- 
cantly lower than the value drawn from the distribution 
above. To account for this effect, the sensitivity predic- 
tion considers the mismatch between the parameters of 
the simulated signals (determined by a random number 
generator) and the template grid points of the search 
(fixed as described earlier). For each simulated source, 
the search grid point that is nearest in the sense of the 
metric is located. Then, using the parameter-space met- 
ric, the mismatch between the simulated signal and the 
closest search template is computed. This gives the frac- 
tional amount by which the 2J^ value is reduced. 

From this ensemble of 2J^ values, one can determine 
the number of coincidences that would be produced by 
each simulated source. As previously described, the post- 
processing sets an effective lower threshold on the J^- 
statistic of the retained candidate events. For each sim- 
ulated source, these thresholds are determined by exam- 
ining the exact workunits that would have contained the 
corresponding signal. Then the number of data segments 
for which the estimated 2J^ values are above threshold is 
counted. If the 2J^ values are above threshold in 12 or 
more of the 17 data segments, the simulated source is 
labeled as "detected", else it is labeled as "undetected". 



D. Search sensitivity, estimated errors, and 
comparison Avith the expected sensitivity 

Shown in Figure [9] are the resulting search sensitiv- 
ity curves as functions of frequency. Each data point on 
the plot denotes the results of 1000 independent trials. 
These show the values of ho as defined in [19 such that 
10%, 50%, and 90% of simulated sources produce 12 or 
more coincidences in the post-processing pipeline. The 
dominant sources of error in this sensitivity curve are 
uncertainties in the magnitudes of the LIGO detector 
response functions (calibration errors). Details of these 
frequency-dependent uncertainties may be found in ref- 
erence [41 . The uncertainties are typically of order 5% 
(LI) and 8% (HI) in the frequency band from 50-1500 
Hz, and are always less than 10%. Systematic errors, 
which arise because of the finite number of Monte-Carlo 



trials and similar effects, are less than ±2%. These can 
be added in quadrature to the uncertainties given in [41 
to obtain frequency-dependent error bounds in the sen- 
sitivity curve. The resulting error in this sensitivity plot 
is below 10% at all frequencies. 

The behavior of the curves shown in Figure |9]essen- 
tially reflect the instrument noise given in Figure^ One 
may fit the curves obtained in Figure [9] to the strain noise 
power spectral density Sh{f) and then describe the three 
sensitivity curves in Figure [9] by 



h^{f)^Rv 



30 hours 



(21) 



where the pre-factors Rx> for different detection proba- 
bilities levels V = 90%, 50%, and 10% are weh fit below 
300 Hz by i?90% = 31.8, i?50% = 20.1, and i?io% = 12.6, 
above 300 Hz by Rqq% = 33.2, i?5o% = 21.0, and Rio% = 
12.9. 

Some other published CW searches were done at 95% 
detection confidence. For comparison in the next section, 
the sensitivity of this search at that confidence is i?95% = 
36.2 below 300 Hz and Rg^% = 37.9 above 300 Hz. The 
iteration method previously described used A/'^^st ~ ^00. 



E. Comparison with sensitivity of other search and 
upper limit methods 

The methods used here would be expected to yield 
very high confidence if a strong signal were present. It is 
interesting to compare the sensitivity of this detection 
method with the sensitivity of CW upper limits such 
as that of reference [27j. The sensitivity of the high- 
confidence detection method used here is well-described 
The same equation describes the re- 



by Equation 21 



suits of the S2 jF-statistic loudest-event upper limit anal- 
ysis [27], but in that work the 95% detection confidence 
curve has a pre- factor Rg^% =29.5. It is useful to under- 
stand the source of this apparent difference in sensitivity 
(a factor of 37.9/29.5 = 1.28). There are three main 
contributors to this difference. 

The most significant difference between the two anal- 
yses is the spacing of the search grid templates. In 
this search, the templates are significantly farther apart 
(worst-case 50% loss of signal-to-noise ratio, or expected 
2J^) than in [27 , where the worst-case mismatch was 
negligible. This effect of employing different mismatches 
has been studied by running the sensitivity estimation 
pipeline using simulated sources only at the template grid 
points, and reduces Rg3% in Equation (21) by a factor of 
1.17. 

Another difference between the two analyses is the de- 
tection criteria. In this work, detection requires a signal 
to produce 12 or more coincidences between the 17 dif- 
ferent data segments. This corresponds to a false alarm 
probability (in Gaussian noise) of order 10~^^ per 0.5 Hz 
frequency-band. This is different from [27], where sim- 
ulated signals are compared against the loudest candi- 
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FIG. 9: Estimated sensitivity of the Einstein@Home search 
for isolated CW sources in the LIGO S4 data set. The set of 
three curves shows the source strain ampUtudes ho at which 
10% (bottom), 50% (middle) and 90% (top) of simulated 
sources would be confidently detected by the Einstein@Home 
pipeline (appearing in coincidence in 12 or more of the 17 data 
segments). The centers of the circles labeled PO to P9 give 
the strain amplitudes of the S4 hardware injections as listed in 
Table [lV| Based on this curve, one would expect that the sim- 
ulated signals P3, P4 and P8 could be confidently detected, 
and that PO, PI and P5 would not be detected. 



date found (largest 2jr). An equivalent detection crite- 
rion for this work would be to compare the simulated 
signals against the loudest candidates (per 0.5 Hz band). 
These typically had 7 or 8 coincidences, corresponding 
to a Gaussian noise false alarm probability of order 10~^ 
and 10~^, respectively. To estimate the effect on the sen- 
sitivity, the sensitivity estimation pipeline was rerun, but 
now requiring the signal to exceed the 2jr-thresholds in 
only 7 of the 17 data segments. This reduced Rq^% in 
Equation (21) by an additional factor of 1.14. 



The least important difference between the two anal- 
yses is the effective threshold on the ^-statistic. As ex- 
plained in Sections |IVA| and |IVC[ only a subset of can- 
didate events with the largest values of 2jr are retained 
in the post-processing, fixing the false alarm probability. 
The smallest 2^- value on this list is typically around 28 
or slightly higher. In [2T a fixed threshold of 2^ = 20 
has been used. Then, events with a combined significance 
below S = 64.4 [see Equation (17)] were also dropped. 
While it is difficult to compare these two criteria, they 
seem to be fairly close. 

Taken together, the differences in grid spacing and de- 
tection thresholds are responsible for, and explain, the 
sensitivity difference in the two analyses (a factor of 
1.17 X 1.14 = 1.33^ 1.28). 



VI. VETOING OF INSTRUMENTAL LINE 
ARTIFACTS 

When the instrument data was prepared and cleaned, 
narrow-band instrumental line features of known origin 
were removed, as previously described in SectionjTTj How- 
ever, the data also contained stationary instrumental line 
features that were not understood, or were poorly under- 
stood. These features were not removed from the input 
data for the search. As a consequence, the output from 
the post-processing pipeline contains instrumental arti- 
facts that in some respects mimic CW signals. But these 
artifacts tend to cluster in certain regions of parameter 
space, and in many cases they can be automatically iden- 
tified and vetoed. In previous incoherent searches for CW 
sources in LIGO data [32] the S-veto method has been 
employed, which excludes the regions of parameter space 
where there is little or no frequency modulation from 
the Earth's motion, leading to a relatively stationary de- 
tected frequency. This cannot directly be applied to a 
coherent matched-filtering search using the ^-statistic. 
Thus the method used here will be similar, but arises 
from a conceptually different approach that is appropri- 
ate for an ^-statistic search. 



A. Parameter space locations of instrumental lines 

For a coherent observation of 30 hours the parameter- 
space regions where instrumental lines tend to appear are 
determined by the global-correlation hypersurfaces [39] in 
parameter space. The global-correlation hypersurface on 
which stationary instrumental lines preferably produce 
candidate events is shown in [49] to be described by 

(U? X Vav) 



/ 



-/(^fiducial) = 



(22) 



where c denotes the speed of light, n is a unit vector 
pointing to the source's sky-location in the SSB frame 
and relates to the equatorial coordinates a and 5 by 
n = (cos J cos a, cos (5 sin a, sin (^), uj is the angular ve- 
locity vector of the Earth as it orbits around the Sun 
(\u)\ ^ 27r/year) and Vav is the average velocity of the 



Earth (Iva 



9.9 X 10- 



This equation can also be 



understood on simple physical grounds. The l.h.s. of 
Equation (22) is the rate of change of detector output 
frequency, for a source whose SSB frequency and spin- 
down are / and /. An instrumental line, which has fixed 
detector frequency, mimics such a source when the l.h.s. 
vanishes. 

The potential CW sources whose locations in parame- 
ter space are consistent with Equation (22) will not pro- 



duce a modulation pattern that would distinguish them 
from an instrumental line. As the resolution in parame- 
ter space is finite, the post-processing analysis eliminates 
(vetoes) candidates that satisfy the condition 
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FIG. 10: The average fraction of sky excluded by the veto 
method as a function of frequency, uniformly averaged over 
the searched spin-down range. 



where the parameter e > accounts for a certain toler- 
ance needed due to the parameter-space gridding. This 
tolerance-parameter can be understood as 



AT 



N, 



cell 



(24) 



where A/ denotes width in frequency (corresponding to 
the coincidence-cell width in the post-processing) up to 
which candidate events can be resolved during the char- 
acteristic length of time AT, and A^ceii represents the size 
of the vetoed or rejected region, measured in coincidence 
cells . In this analysis AT = 2122 735 s {^ 24 days) is 
the total time interval spanned by the data 



AT 
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end 



min tn 



r 



end 
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(25) 



For potential sources that satisfy (23), the modulation 



due to the Earth's motion does not make the signal ap- 
pear in more than Nceii coincidence cells during AT. 



veto depends upon whether or not (as the source po- 
sition varies over the sky) the excluded region between 
the lines lies inside or outside of the boundaries of the 
search, or intersects it. In this work, for the search re- 
gion —f/r < f < 0.1 //r described in the Abstract, 
the excluded region lies entirely within the parameter 
space above 300 Hz, but crosses the boundaries below 
300 Hz. This is because a wider range of spin-down ages 
is searched below 300 Hz. 

The fractional volume of the region in parameter space 
excluded by the veto may be easily calculated. The de- 
tails of the calculation are in Appendix [Aj The resulting 
fraction of sky excluded by the veto (uniformly averaged 
over spin-down) as a function of frequency is shown in 
Figure [lO] In this search, the fraction of the sky ex- 
cluded for frequencies / G [300, 1500) Hz has been fixed 
at the constant fraction 30%. In this search, the frac- 
tion of the sky excluded for frequencies / G [50, 300) Hz 
has been chosen to depend upon the values of / and /, 
where the uniform average of the excluded sky fraction 
over the spin-down range considered in this analysis is 
36% at 50 Hz and 6% just below 300 Hz. Finally, Fig- 
ure [15] shows a conclusion diagram illustrating which of 
the candidates have been vetoed in this search. 



VII. HARDWARE-INJECTED SIGNALS 

A good way to test and validate search algorithms and 
code is to add simulated signals into the detector's data 
stream. This can either be done while the experiment 
is in progress (real-time injections) or after the data has 
been collected (software injections). If it is done while 
the experiment is in progress, the simulated signals can 
either be added into the hardware (into feedback and 
error-point control signals) or after data acquisition. 

At the time that the S4 run was carried out, ten sim- 
ulated CW signals were injected at the hardware level: 
using magnetic coil actuators, the interferometer mirrors 
were physically made to move as if a gravitational wave 
was present. 



A. Parameters of hardAvare injections 



B. Fraction of parameter space excluded by the 
veto 

One can visualize and calculate the volume of the re- 
gion in four-dimensional parameter space that is excluded 
by this veto. For a given source sky position. Equa- 
tion (23) is linear in / and /. Thus, for fixed sky position 



n, the veto condition (23) defines two parallel lines in the 



(/, /)-plane. Candidate events that lie in the region be- 
tween the lines are discarded (vetoed). Candidates that 
lie outside this region are retained (not vetoed). The lo- 
cations of these two lines in the (/, /) plane depend upon 
the sky position. The fractional volume excluded by the 



Table IV shows the parameters of the hardware injec- 
tions that were carried out at the LIGO detectors during 
the S4 run, mimicking gravitational-wave signals from 
ten different isolated pulsars with different frequencies, 
sky locations, and frequency derivatives. The ten ar- 
tificial pulsars are denoted PulsarO to Pulsar9. At the 
time of the injections, lack of complete knowledge of the 
instrument's response function (calibration) meant that 
the actual hardware injections did not actually have the 
intended strain amplitudes as given in the Table. The ef- 
fective strain amplitudes may be computed from correc- 
tion factors provided in reference [32]. These factors are 
1.12 for all simulated pulsars in the HI detector. In the 



18 




^^^^^^. 


a 


^ 




^ 


a 


n^ 




> 


a 


H 


ii 


f^ 


a 


^^^^^^ 




^ 
B 




FIG. 11: Einstein@Home results showing the 10 hardware- injected pulsar signals labeled PO to P9. Here, a narrow band of 
width 2 X lO"'^ / to either side of each injection's frequency / is considered. The color-bar in each plot indicates the number of 
coincidences. As shown in the color-scale, only candidates having 7 or more coincidences appear. For each hardware injection 
a group of three different sub-plots are given representing different projections of the parameter space. The left sub plot is a 
Hammer- Aitoff projection of the entire sky. The middle sub-plot shows declination S versus frequency /. The right sub-plot 
shows spin-down / versus frequency /, where the region between the two horizontal magenta dashed lines refers to searched 
range of spin-downs. The center of a magenta circle represents the location of the injection. P4 and P8 appear at the wrong sky 
position because their intrinsic spin-downs lie outside the searched range. Table [Vl] shows a comparison with the expectations 
for these simulated signals. 



LI detector, the correction factor is 1.11 for all simulated 
pulsars, except for Pulsar 1 (1.15) and Pulsar9 (1.18). 

B. Duty cycle of hardAvare injections 

During S4 the hardware injections were not active all 
of the time. Table [V| shows the fractional overlap between 
the times when the hardware injections were active and 
the times of the S4 Einstein@home data segments. As 
can be seen from the table, the hardware injections were 
only turned on during twelve of the data segments ana- 
lyzed in this paper, and for two of those twelve data data 
segments, the injections were only turned on for about 
20% of the data taking time. In the remaining ten data 
segments, the hardware injections were turned on for al- 



most the entire segment. This needs to be taken into 
account when analyzing the Einstein@Home search re- 
sults for these injections. Because of this, the maximum 
possible number of coincidences expected from these sim- 
ulated signals is 12, even though 17 data segments are 
searched. 



C. Results from the hardAvare injections 



For each hardware-injected pulsar signal Table [VI| com- 
pares a prediction for the outcome of the Einstein@Home 
search to the actual results found through the Ein- 
stein@Home analysis pipeline. The predicted values 
given in Table [Vl] are obtained by feeding the sensitivity- 
estimation pipeline, which was described in Section |V| 
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Name 


/(tref) [Hz] 


/ [Hz s" 


-'] 


a [rad] 


6 [rad] 


ip [rad] 


$0 [rad] 


cos^ [rad] 


ho 




PulsarO 


265.57693318 


-4.15 X 10" 


-12 


1.248817 


-0.981180 


0.770087 


2.66 


0.794905 


4.93 X 10" 


-25 


Pulsarl 


849.07086108 


-3.00 X 10" 


-10 


0.652646 


-0.514042 


0.356036 


1.28 


0.463799 


4.24 X 10" 


-24 


Pulsar2 


575.16356732 


-1.37 X 10" 


-13 


3.756929 


0.060109 


-0.221788 


4.03 


-0.928575 


8.04 X 10" 


-24 


Pulsars 


108.85715940 


-1.46 X 10" 


-17 


3.113189 


-0.583579 


0.444280 


5.53 


-0.080666 


3.26 X 10" 


-23 


Pulsar4 


1402.11049084 


-2.54 X 10" 


-08 


4.886707 


-0.217584 


-0.647939 


4.83 


0.277321 


4.56 X 10" 


-22 


PulsarS 


52.80832436 


-4.03 X 10" 


-18 


5.281831 


-1.463269 


-0.363953 


2.23 


0.462937 


9.70 X 10" 


-24 


PulsarG 


148.44006451 


-6.73 X 10" 


-09 


6.261385 


-1.141840 


0.470985 


0.97 


-0.153727 


2.77 X 10" 


-24 


Pulsar? 


1220.93315655 


-1.12 X 10" 


-09 


3.899513 


-0.356931 


0.512323 


5.25 


0.756814 


1.32 X 10" 


-23 


Pulsars 


193.94977254 


-8.65 X 10" 


-09 


6.132905 


-0.583263 


0.170471 


5.89 


0.073904 


3.18 X 10" 


-23 


Pulsar9 


763.8473216499 


-1.45 X 10" 


-17 


3.471208 


1.321033 


-0.008560 


1.01 


-0.619187 


8.13 X 10" 


-24 



TABLE IV: Parameters for hardware-injected CW signals during the S4 run, labeled PulsarO to Pulsar9. The parameters 
are defined at the GPS reference time tref = 793130413s in the Solar System Barycenter. These are the frequency /(tref), the 
spin-down /, the sky position right ascension a and declination S, the polarization angle ip, the initial phase $o, the inclination 
parameter cos^, and the dimensionless strain amplitude ho. Because the calibration was only accurately determined after 
S4 was finished, the HI strain amplitudes should be multiplied by the correction factor 1.12. The LI amplitudes should be 
multiplied by 1.15 for Pulsarl, 1.18 for Pulsar9, and 1.11 for the others. 



j 


Detector 


Overlapping 


Fractional 






Duration [s] 


Overlap 


1 


HI 


107201 


99.3 % 


2 


HI 


107554 


99.6 % 


3 


HI 


107272 


99.3 % 


4 


HI 








5 


HI 


99799 


92.4 % 


6 


HI 








7 


HI 


101991 


94.4 % 


8 


HI 


21268 


19.7 % 


9 


HI 


100773 


93.3 % 


10 


HI 








11 


LI 


23164 


21.5 % 


12 


LI 


106760 


98.9 % 


13 


LI 


107294 


99.4 % 


14 


LI 


102711 


95.1 % 


15 


LI 








16 


LI 








17 


LI 


98696 


91.4 % 



TABLE V: The time overlap between the Einstein@Home 
data segments and the hardware injections. The hardware 
injections were only turned on about 2/3 of the time. The 
columns are data segment index j, detector, the duration of 
the overlap, and the fractional overlap (obtained by dividing 
the third column by 30 hours = 108 000 s). 



with the parameters of the simulated pulsars and only 
considering data segments where the hardware injections 
were active. 

As shown in Table [Vl] and consistent with Figure [9) the 
hardware-injected signals PulsarO, Pulsarl, PulsarS and 
PulsarG are too weak to be confidently detected by the 
search. In contrast, Pulsar2, Pulsar3, Pulsar4 and Pul- 



Name 


Predicted Measured Predicted Measured 




C 


C 


S 


S 


Pulsar2 


12 


13 


263.1 


249.3 


Pulsar3 


12 


12 


3160.9 


2397.5 


Pulsar 4 


12 


12 


35108.2 


1749.6 


Pulsar7 


6 


7 


93.2 


100.0 


Pulsar8 


12 


13 


3692.6 


2263.6 


Pulsar9 


7 


7 


131.2 


98.9 



TABLE VL The estimated (predicted) and obtained (mea- 
sured) results for the hardware-injected pulsar signals. For 
each simulated signal the predicted number of coincidences C 
and a predicted value for the significance S is given, as well as 
the measured number of coincidences and measured value for 
the significance from the Einstein@Home search. The mea- 
sured values are obtained by maximizing over a narrow band 
of 2 X 10~^ / on either side of the injection frequency, the 
whole sky and the entire spin-down range. As explained in 
the text, Pulsar4 and Pulsar8 are not expected to have the 
correct significance. PulsarO, Pulsarl, Pulsar5 and Pulsar6 
are not listed. They are so weak that they produce less than 
7 coincidences, consistent with random noise containing no 
signal. 



sarS are clearly detected. The parameters of Pulsar7 and 
PulsarQ are such that in both cases the search pipeline 
found 7 coincidences, but this is consistent with the level 
of coincidences that would result from Gaussian noise 
with no signal present, and so these are not confidently 
detected. 

Figure pT] presents the results of the search for all hard- 
ware injections. Small subspaces of the search parame- 
ter space around the hardware injections are shown, as 
well as the locations of the artificial signal parameters. 



The subspaces considered in Figure 11 and also for the 
(measured) results presented in Table VI are constrained 
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FIG. 12: Numbers of coincidences greater than 7 (top) and 
the significance (bottom) of all candidates found in the Ein- 
stein@Home post-processing, shown as functions of frequency. 
The light-gray shaded rectangular regions highlight the S4 
hardware injections, listed in Table |IV| The data points col- 
ored in dark-gray show the candidates resulting from the 
hardware- injected CW signals. 



to a band of 2 x 10~^ / to either side of the injected 
frequency. This choice of frequency-bandwidth is moti- 
vated by the global parameter-space correlations [49l [50] 
and represents the approximate maximum extension in 
frequency-direction of a global-correlation hypersurface 
generated by a single signal. 

Moreover, the global correlations also explain the sig- 
nificant sky position offset between the a priori location of 
the simulated source and the location where the search 
located the source with respect to the detected signals 
Pulsar4 and PulsarS. This arises because for Pulsar4 and 
PulsarS, the spin-down range that is searched (region be- 
tween dashed lines in the far right column) is too small 
to include the actual spin-down value used in creating 
the simulated signals. Due to the global parameter-space 
correlations the offset between the actual and detected 
spin-down value gives rise to the offset in the sky posi- 
tion. The observed structure of large-coincident events 
in the sky is consistent with a (sky-projected) global- 
correlation hypersurface first found analytically in [49] . 
This is also why Pulsar4 shows a considerable discrep- 
ancy between the significance that would have been ex- 
pected if the search-grid had also covered the a priori 
parameters, and the significance that was actually ob- 



served in the search, as shown in Table VI 



VIII. RESULTS 



This section presents the results of the Einstein@Home 
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FIG. 13: All candidates obtained from the post-processing 
that have more than 7 coincidences, shown in Hammer- Aitoff 
projections of the sky. The color-bar indicates the number 
of coincidences of a particular candidate (cell). The upper 
plot includes the S4 hardware-injected pulsars. In the lower 
plot, bands of 2 x 10~^ / width to either side of the hardware 
injections' frequencies / have been removed. 



post-processing results, from 50 to 1500 Hz. In Figure 12 
the coincidences and significance of all candidates that 
have 7 or more coincidences are shown as functions of 
frequency. Figure [13| presents the same information as 
given in Figure [12] but projected on the sky, and showing 
all cells that have more than 7 candidate events. 



In Figure 13 the number of coincidences is maximized 
over the entire sky and full spin-down range. The color 
indicates the numbers of coincidences, where the same 
color-scale has been used in each plot. The maximum 
possible number of coincidences ranges from a minimum 
of to a maximum of 17 (the number of data segments 
analyzed). The meaning of coincidences is that there 
is no candidate event found, 1 coincidence means a sin- 
gle candidate events is found (which is always coincident 
with itself). 

Four illustrative examples of different types of typi- 
cal post-processing results in particular 10 Hz bands are 
shown in Appendix [B] by Figures [T6| [18] [17] and [19] 



S4 CW search. Figures 12 and 13 give a summary of all 



Table VII shows all candidates (cells) which have 10 
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/cand [Hz] /start [Hz] A/eand [Hz] 


(^cand [rad] 


QJcand [rad] 


/cand [Hz S" 


-'] 


'-'cand 


^Hl 

<-cand 


r^Ll 
^cand 


*Scand 


Information 


193.9276 


193.9263 


0.040112 


-0.583514 


4.723595 


-5.6001 X 10- 


-09 


13 


7 


6 


2263.6 


Pulsar 8 


575.1681 


575.1562 


0.030612 


0.285505 


3.834511 


-5.0913 X 10" 


-10 


13 


7 


6 


249.3 


Pulsar 2 


1128.1147 1128.0336 


0.220321 


-1.395918 


0.744273 


-3.4249 X 10" 


-09 


13 


10 


3 


219.3 


HI MC 2/3 


108.8549 


108.8522 


0.008158 


-0.705729 


3.361465 


-4.4362 X 10" 


-11 


12 


9 


3 


2397.5 


Pulsar 3 


329.6107 


329.5507 


0.066447 


1.027320 


1.336051 


-5.7799 X 10" 


-10 


12 


10 


2 


3127.1 


Demod 


545.9973 


545.9929 


0.10958 


-0.293877 


4.849960 


-1.5782 X 10" 


-09 


12 


10 


2 


893.3 


HI MC 2/3 


566.0868 


566.0490 


0.105853 


-1.367663 


0.665233 


-1.626 X 10" 


-09 


12 


10 


2 


2340.8 


HI MC 2/3 


568.0886 


567.9893 


0.165769 


-1.323532 


0.726729 


-1.7149 X 10" 


-09 


12 


10 


2 


4137.7 


HI MC 2/3 


648.8288 


648.6930 


0.206223 


-1.232868 


1.005733 


-1.0298 X 10" 


-09 


12 


10 


2 


1870.8 


HI MC 1 


1143.9976 


1143.9182 


0.232221 


-1.491264 


1.314456 


-7.7434 X 10" 


-10 


12 


10 


2 


1028.8 


HI Cal 


1144.5198 


1144.4533 


0.228407 


-1.535248 


4.497733 


-2.5257 X 10" 


-11 


12 


10 


2 


989.8 


HI Cal 


1289.6769 


1289.5081 


0.242915 


1.461093 


0.266878 


-2.0949 X 10" 


-09 


12 


10 


2 


493.7 


HI MC 1 


1402.2838 


1402.2677 


0.063117 


1.025583 


2.502838 


-3.8482 X 10" 


-09 


12 


6 


6 


1749.6 


Pulsar 4 


329.7593 


329.7396 


0.066078 


-1.536179 


4.887048 


-5.5375 X 10" 


-10 


11 


10 


1 


3038.3 


Demod 


335.7735 


335.7100 


0.065415 


0.469606 


0.955884 


-1.0646 X 10" 


-09 


11 


10 


1 


298.5 


EM Interference 


545.9232 


545.8662 


0.063608 


-1.060735 


1.078303 


-8.032 X 10" 


-10 


11 


10 


1 


196.8 


HI MC 2/3 


564.1219 


564.0096 


0.113783 


0.386877 


1.111355 


-1.6868 X 10" 


-09 


11 


10 





1069.3 


HI MC 2/3 


646.3758 


646.3206 


0.127884 


-1.281366 


0.897933 


-1.8931 X 10" 


-09 


11 


10 


1 


3202.7 


HI MC 1 


1092.1387 1091.9671 


0.217482 


-0.523866 


1.302500 


-6.4347 X 10" 


-11 


11 


10 


1 


196.7 


HI MC 2/3 


1136.2217 1136.1460 


0.168345 


-1.216945 


0.935876 


-3.4811 X 10" 


-09 


11 


10 


1 


165.6 


HI MC 2/3 


1142.8210 


1142.7200 


0.23173 


-1.310037 


1.114563 


-3.5022 X 10" 


-09 


11 


10 


1 


250.7 


HI Cal 


1145.8318 


1145.6515 


0.231067 


1.330065 


0.976422 


-2.0297 X 10" 


-10 


11 


10 


1 


256.4 


HI Cal 


1376.7370 


1376.4697 


0.271536 


0.201677 


1.282354 


-2.3875 X 10" 


-09 


11 


9 


2 


165.0 


TM violin 


1388.6402 


1388.4070 


0.279967 


1.176082 


0.850794 


-2.8907 X 10" 


-09 


11 


10 


1 


200.0 


TM violin 


56.9966 


56.9966 




-0.935903 


0.150238 


-1.5029 X 10" 


-09 


10 


8 


2 


136.7 


EM Interference 


329.4918 


329.4784 


0.021358 


-1.307440 


4.692056 


-5.2405 X 10" 


-10 


10 


10 





1137.8 


Demod 


392.8322 


392.8322 




-1.210088 


1.268596 


-1.069 X 10" 


-09 


10 


9 


1 


150.9 


HlCal 


393.4060 


393.4057 


0.000342 


0.632053 


1.270922 


-1.1043 X 10" 


-09 


10 


9 


1 


154.7 


HlCal 


646.7224 


646.7224 


0.002174 


-1.446520 


0.825633 


-1.8813 X 10" 


-09 


10 


3 


7 


2774.4 


LI MC 1 


648.4132 


648.4132 


0.024291 


1.319729 


1.033730 


-1.8479 X 10" 


-09 


10 


3 


7 


5067.5 


LI MC 1 


658.6353 


658.6353 


0.000055 


-0.470832 


4.762475 


-1.6992 X 10" 


-09 


10 


3 


7 


261.4 


EX +15v 


777.9202 


777.8377 


0.117087 


1.511859 


4.010213 


-5.6101 X 10" 


-10 


10 


3 


7 


1951.7 


EM Interference 


1296.4962 


1296.4962 




-0.993190 


4.557370 


-1.0022 X 10" 


-09 


10 


3 


7 


247.1 


LI MC 1 



TABLE VII: The post-processing candidates that have 10 or more coincidences. The frequency /cand corresponds to the most 
coincident candidate in the band. The lowest frequency of a candidate in the band is labeled by /start- The difference from 
the highest frequency is given by A/cand- The parameters ^cand, Q^cand, /cand, Ccand = Cf^nd + ^Jknd and <Scand are for the most 
significant most coincident candidate within the frequency band, where C^^nd and Cj'and denote the number of coincidences 
contributing to Ccand from detector HI and LI, respectively. The column "Information" lists information about the source. 
The following are understood sources of narrow-band line noise in the instrument: "Demod" are the electronics boards that 
demodulate the signal at the antisymmetric port of the interferometer, "HI (or LI) MC 1" is a violin mode resonance of the 
first mode cleaner mirror, "HI MC 2/3" are violin mode resonances of the second and third mode cleaner mirrors, "TM violin" 
are harmonics of the test mass violin modes, "EX H-15v" is a fifteen volt power supply at the end station of the X arm, "EM 
Interference" is electromagnetic interference, "HI Cal" are side-bands of calibration lines at 393.1Hz and 1144.3 Hz. 



or more coincidences. In cases where a set of candidates 
is clustered together at slightly different frequencies, Ta- 
ble |VTT] lists the bandwidth in frequency covered by these 
candidates and shows the parameters of the most co- 
incident candidate. If candidates within these narrow 
frequency-bands have the same number of coincidences, 
then the candidate with the largest significance is shown. 



Table VIII shows the same information after the veto 



method described in Section [VI] has been applied, for 
candidates with 9 or more coincidences. There are no 
candidates that exceed the predefined detection thresh- 
old of appearing in 12 or more data segments. (Note that 
this would be a threshold for initiating a more extensive 
investigation of the candidate event, not a threshold for 
announcing a discovery!) 

Figure [l4] shows all candidates from the post- 
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/cand [Hz] /start [Hz] A/eand [Hz] 


Scand [rad] 


«cand [rad] 


/cand [Hz S-1] 


^cand 


'-'cand 


'-'cand 


*Scand 


Information 


193.9276 


193.9261 


0.040646 


-0.583514 


4.723595 


-5.6001 X 10"°^ 


13 


7 


6 


2263.6 


Pulsar 8 


575.1681 


575.1562 


0.039394 


0.285505 


3.834511 


-5.0913 X 10-^° 


13 


7 


6 


249.3 


Pulsar 2 


108.8549 


108.8518 


0.008506 


-0.705729 


3.361465 


-4.4362 X 10"^^ 


12 


9 


3 


2397.5 


Pulsar 3 


1402.2838 


1402.2488 


0.08678 


1.025583 


2.502838 


-3.8482 X 10"°^ 


12 


6 


6 


1749.6 


Pulsar 4 


545.9987 


545.9568 


0.141563 


-0.398855 


5.013332 


-4.6693 X 10"^° 


11 


10 


1 


794.1 


HI MC 2/3 


56.9966 


56.9963 


0.000933 


-0.935903 


0.150238 


-1.5029 X 10-°^ 


10 


8 


2 


136.7 


EM Interference 


329.4849 


329.4833 


0.005843 


-0.344739 


5.171401 


-5.6694 X 10-^° 


10 


10 





1024.0 


EM Interference 


329.6040 


329.6040 




-0.439100 


1.006331 


-5.2546 X 10"^° 


10 


9 


1 


2625.6 


EM Interference 


329.7434 


329.7413 


0.032463 


-0.338712 


5.025108 


-5.6923 X 10"^° 


10 


9 


1 


2490.4 


EM Interference 


567.9984 


567.9984 


0.051768 


-0.353846 


5.116972 


-1.5532 X 10"°^ 


10 


9 


1 


409.3 


EM Interference 


69.6964 


69.6964 




-1.223613 


4.232687 


-5.4823 X 10"^° 


9 


9 





130.3 


EM Interference 


317.4207 


317.4207 




1.389330 


2.663214 


-8.0338 X 10-^° 


9 


3 


6 


157.8 


EM Interference 


329.5615 


329.5615 




-1.027976 


3.822726 


-6.3014 X 10"^° 


9 


7 


2 


2176.0 


Demod 


335.7541 


335.7141 


0.056927 


1.395059 


3.271989 


-6.362 X 10"^° 


9 


9 





259.3 


EM Interference 


795.4783 


795.4783 




0.245291 


3.211417 


-1.4374 X 10"°^ 


9 


7 


2 


110.7 
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1117.3032 
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116.1 


Unknown 


1145.6678 
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9 
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168.7 
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TABLE VIII: Post-processing candidates that have 9 or more coincidences and that are not excluded by the veto. The frequency 
/cand corresponds to the most coincident candidate in the band. The lowest frequency of a candidate in the band is labeled by 
/start. The difference from the highest frequency is given by A/cand- The parameters ^cand, Q^cand, /cand, Ccand = CfJ^d + ^^J'and 
and tScand are for the most significant most coincident candidate within the frequency band, where C^^nd ^^^ ^Jknd denote 
the number of coincidences contributing to Ccand from detectors HI and LI, respectively. The column "Information" lists 
information about the source. The following are understood sources of narrow-band line noise in the instrument: "Demod" are 
the electronics boards which demodulate the signal at the antisymmetric port of the interferometer, "HI MC 2/3" are violin 
mode resonances of the second and third mode cleaner mirrors, "EM Interference" is electromagnetic interference, "HI Cal" 
are side-bands of a 1144.3 Hz calibration line. For the single candidate labeled "Unknown" in the last column no instrumental 
source could be confidently identified, however the 9 coincidences are far below the confident-detection threshold. 



processing results that have not been discriminated by 



the veto introduced in Section VI Figure 15 illustrates 
the fraction of candidates that has been excluded by the 
veto. Removing fractional bands of2xl0~^ / around the 
frequencies / of the S4 hardware injections, the veto dis- 
criminates 99.5% of all candidates that have more than 
7 coincidences. 



IX. CONCLUSION 

These are the first published results from the Ein- 
stein@Home project, which was launched in February 
2005. While no credible CW sources were found in this 
search of LIGO S4 data, the results clearly establish that 
this type of distributed computing project can carry out 
a credible and sensitive search for such signals. 

In retrospect, it probably would have been a good idea 
to employ identical grids on the four-dimensional param- 
eter space for all 17 data segments. This would have 
required more CPU time on the part of participants, but 
would have greatly simplified and sped up the develop- 
ment of the post-processing pipeline and would also have 
greatly simplified the interpretation of the results. 

A similar search (also with a 30-hour time baseline) 



has already been completed using 660 hours of data from 
the beginning of the S5 science run. The post-processing 
of that data set is currently underway, using methods 
identical to those employed here. 

Future Einstein@Home searches overcome some of the 
sensitivity limitations discussed at the end of Section |V| 
by doing the incoherent step (called "post-processing" in 
this paper) on the host machines. This allows the use 
of the optimal threshold of 2J^ ~ 5, so those searches 
are expected to be the most sensitive blind CW searches 
that will be possible using LIGO data. Results from those 
searches should become available within the next one to 
two years, and are expected to offer more than one order 
of magnitude improvement in strain sensitivity compared 
with the work presented here. 

In the longer term, further increases in sensitivity will 
probably result mostly from improvements in the detec- 
tors rather than from improvements in the data anal- 
ysis methods. In 2009 LIGO is expected to begin its 
S6 run with an "enhanced" detector configuration that 
should improve on S5 sensitivity by at least a factor of 
two. By 2014, an advanced LIGO detector configuration 
should give at least another factor of five improvement. 
By combining these data sets with those from LIGO's 
international partner projects Virgo and GEO, there is 
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FIG. 14: Candidates not eliminated by the veto. This 
shows Hammer- Ait off sky projections of all candidates ob- 
tained from post-processing that had more than 7 coinci- 
dences and that passed the veto. The upper plot includes 
the S4 hardware injections. The lower plot removes bands 
of 2 X 10~^ / width to either side of the hardware injections' 
frequencies /. In comparison to Figure [13] after excluding the 
hardware injections, the veto rejects 99.5% of all candidates. 



real hope that the first direct CW detection can be made 
using methods like the ones described here. 
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FIG. 15: Conclusion diagram of candidates discriminated 
by the veto method. All candidate cells obtained from post- 
processing that have more than 7 coincidences are shown, 
where the color-bar indicates the number of coincidences of a 
particular cell. The vertical axis represents the veto quantity 
on the left-hand side of (23), as a function of frequency. Can- 



didates located below the magenta line are eliminated by the 
veto. The four accumulations of highly coincident cells above 
the magenta line are the hardware injected pulsars, which are 
not eliminated by the veto. 
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The fractional volume of the region in parameter space 
excluded by the veto method presented in [49 and used in 
Section [Vl| may be easily calculated. Since the time AT 
is small compared to one year, one may use the following 
approximation 



{iV X Va 



(Al) 



where G [0, tt] is the angle between the SSB-to-Earth 
vector and the source sky position n. The veto condi- 
tion (23) may then be rewritten as 



\f + jfcosO\<e 



(A2) 
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where 7 is defined as 7 = |cj||vav|/c. For fixed values of 
/ and / the situation is depicted in Table IX Depending 
upon the values of (±e — /)/7/, a part of the sky might 
be excluded by the veto. As shown in the Table, there 
are six possible cases, determined by the values of 



cos^_ 
cos6>+ 



7/ 
7/ " 



and 



(A3) 

(A4) 



For example in the case (labeled case 4 in Table IX ) where 
both cos^_ and cos6>+ lie in the range [— l,l]then the 
excluded region of the sky is an annulus defined by < 
0- < < 0-^ < TT^ and the excluded solid angle is 



a 



excluded 



27r 



d(j) 



dO^. 



= 27r(cos6>_ — cos^+) . 
The fraction of sky excluded in this case is then 

^ ^excluded 



47r 



= (cos i9_ -cosi9+)/2 



7/ 



(A5) 
(A6) 

(A7) 
(A8) 

(A9) 



In the other cases listed in Table [IX| the excluded region 
of the sky might be a cap about 6> = or about 6> = tt or 
the null set, or the entire sky. 

In this search, the fraction of the sky excluded for fre- 
quencies / G [300, 1500) Hz has been fixed at the con- 
stant fraction ^excluded /47r = 30%. This is equivalent to 
choosing e to be a linear function of frequency 



1 

7 



0.3^ 



/ 



(AlO) 



In this search, the fraction of the sky excluded for fre- 
quencies / G [50, 300) Hz has been chosen to depend upon 
the value of /. The instruments allow (e.g. compare 
with Figure fT5| the use of a frequency-independent value 
e = 5.4 X 10^^^ Hz/s which corresponds to A^ceii = 1-5. 
Within the region of parameter space which is searched 
(— //r < / < 0.1 //r for r = 1000 years) cases 4, 5, or 
6 from Table [IX| occur depending of the spin-down value 

/. If 



/>e 



l<^l|Vav|/ 



(All) 



then case 4 of Table |IX| applies, and the fraction of sky 
excluded is given by 



Vt, 



excluded 



4^ -^|u.||Vav|/- ^^^^^ 

This fraction ranges from 52% at 50 Hz to 8.7% at 300 Hz. 
If / is in the interval 



/e 



|a;|Kv|/ lu'lKvl/ 



(A13) 



then case 5 of Table |IX] applies, and the fraction of sky 
excluded is given by 



Q.^ 



excluded 

47r 



1 



/ 



Finally, if 



/< 



/|w||Vav|/c 



l<^l|Vav|/ 



(A14) 



(A15) 



then case 6 applies and none of the sky is excluded 
by the veto: ^excluded = 0. Below 300 Hz, one can 
compute a uniform average of the excluded sky fraction 
over the spin-down range considered in this analysis. As 



shown in Figure 10 this gives an excluded sky fraction 
of 36% at 50 Hz and 6% just below 300 Hz. 



APPENDIX B: ILLUSTRATIVE EXAMPLES OF 
TYPICAL POST-PROCESSING RESULTS 



In the following, illustrative examples of different types 
of typical post-processing results in four individual 10 Hz 
bands are shown in Figures [16) [18) [17) and 19 Fig- 
ure 



16 shows a 10 Hz frequency band containing pure 



Gaussian noise. Figure [18] shows the frequency band of 
the hardware- injected signal Pulsar2. Figure [T7| shows a 
"quiet" 10 Hz band of real instrument data without any 
"noisy" lines. In contrast to this. Figure [19] shows a noisy 
band which is polluted by instrumental noise artifacts. 

As described in Section [VIII[ each of these plots shows 
the number of coincidences maximized over the entire 
sky and full spin-down range. The color indicates the 
numbers of coincidences, where the same color-scale has 
been used in each figure. 
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FIG. 16: Einstein@Home S4 Post-Processing results for the 
frequency band 340.0-350.0 Hz, which is pure Gaussian noise 
for LI and mostly Gaussian noise for HI. This is because in 
this band the line-cleaning process has replaced all the LI data 
and most of the HI data with computer-generated random 
numbers (see Table[Il|. From top to bottom the different plots 
show the numbers of coincidences in a 3D map of sky and 
spin-down, in a 2D Hammer-Aitoff projection of the sky, in a 
2D plot of declination over frequency, and in a histogram as 
a function of frequency. 
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FIG. 17: Einstein@Home S4 Post-Processing results for a 
"quiet" frequency band of real instrumental data from 110.0- 
120.0 Hz. From top to bottom the different plots show the 
numbers of coincidences in a 3D map of sky and spin-down, 
in a 2D Hammer- Aitoff projection of the sky, in a 2D plot of 
declination over frequency, and in a histogram as a function 
of frequency. 



FIG. 18: Einstein@Home S4 Post-Processing results for the 
frequency band 570.0-580.0 Hz including a hardware injected 
CW signal (Pulsar2). From top to bottom the different plots 
show the numbers of coincidences in a 3D map of sky and 
spin-down, in a 2D Hammer- Aitoff projection of the sky, in a 
2D plot of declination over frequency, and in a histogram as 
a function of frequency. 
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FIG. 19: Einstein@Home S4 Post-Processing results for a 
"noisy" frequency band of data polluted by instrumental noise 
artifacts from 640.0-650.0 Hz. These spectral features are res- 
onance modes of the mode cleaner optics suspensions. From 
top to bottom the different plots show the numbers of coinci- 
dences in a 3D map of sky and spin-down, in a 2D Hammer- 
Aitoff projection of the sky, in a 2D plot of declination over 
frequency, and in a histogram as a function of frequency. 
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